Superfluid Bose-Fermi mixture from weak-coupling to unitarity 
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We investigate the zero-temperature properties of a superfluid Bose-Fermi mixture by introduc- 
ing a set of coupled Galilei-invariant nonlinear Schrodinger equations valid from weak-coupling to 
unitarity. The Bose dynamics is described by a Gross-Pitaevskii-type equation including beyond- 
mean-field corrections possessing the correct weak-coupling and unitarity limits. The dynamics of 
the two-component Fermi superfluid is described by a density-functional equation including beyond- 
mean-fleld terms with correct weak-coupling and unitarity limits. The present set of equations is 
equivalent to the equations of generalized superfluid hydrodynamics, which take into account also 
surface effects. The equations describe the mixture properly as the Bose-Bose repulsive (positive) 
and Fermi- Fermi attractive (negative) scattering lengths are varied from zero to infinity in the pres- 
ence of a Bose-Fermi interaction. The present model is tested numerically as the Bose-Bose and 
Fermi-Fermi scattering lengths are varied over wide ranges covering the weak-coupling to unitarity 
transition. 

PACS numbers: 03.75.Ss, 03.75.Hh 
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I. INTRODUCTION 

The macroscopic quantum phenomenon of superflu- 
idity has been clearly demonstrated in recent experi- 
ments with ultracold and dilute gases made of alkali- 
metal atoms [l|, Q • Both bosonic and fermionic ultracold 
and dilute superfluids can be accurately described by the 
zero-temperature hydrodynamical equations of superflu- 

ids flaiii. 

Trapped Bose-Fermi mixtures, with Fermi atoms in a 
single hyperfine state (normal Fermi gas), were investi- 



there was a study of a dilute mixture of superfluid bosons 
and fermions across a Feshbach resonance of the Fermi- 
Fermi scattering length a f , obtaining the phase diagram 
of the mixture in a box [1^ . In the strict one-dimensional 
case, we found that the superfluid Bose-Fermi mixture 
exhibits phase separation and solitons by changing the 
Bose-Fermi interaction [l^ [2^ . 

In the present paper we analyze in detail a 3D su- 
perfluid Bose-Fermi mixture under harmonic trapping 
confinement when the Bose-Bose repulsive (positive) and 
Fermi-Fermi attractive (negative) scattering lengths are 
varied from very small (weak coupling) to infinitely large 
(strong coupling) values. The variation of these scatter- 
ing lengths from weak to strong coupling is achieved in 
laboratory pl| by varying a background magnetic field 
near a Feshbach resonance. 

In the first part of the paper we derive Galilei-invariant 
nonlinear Schrodinger equations for superfluid Bose and 
Fermi dynamics, valid from weak-coupling to unitarity 
in each case, which are equivalent to the hydrodynam- 
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ical equations. For bosons the equation is a general- 
ized [22,, 23, 24| zero-temperature Gross-Pitaevskii (GP) 
equation [25| including beyond-mean-field corrections to 
incorporate correctly the effect of bosonic interaction for 
large and positive (repulsive) scattering lengths ai,. In 
the ordinary GP equation the nonlinear term is linearly 
proportional to at and hence highly overestimates the 
effect of bosonic interaction for large positive values of 
ab- In fact there is a saturation of bosonic interac- 
tion in the so-called unitarity at —>■ +00 limit, prop- 
erly taken care of in the present generalized GP equa- 
tions. For superfluid fermions the present formulation is 
based on the density- functional (DF) theory [26|, H?! for 
fermion pairs [H, [2^, H^, [13, HH including beyond- mean- 
field corrections to incorporate properly the effect of 
fermionic attraction between spin-up and down fermions 
specially for large negative (attractive) values of Fermi- 
Fermi scattering length a/. The DF theory, in different 
forms, has already been applied to the problem of ultra- 
cold fermions [12, [1^ • The usual DF equation is valid 
for the weak-coupling Bardeen-Cooper-Schreiffer (BCS) 
limit: a/ —0. The present generalized DF equation 
correctly accounts for the dynamics as the Fermi-Fermi 
attraction is varied from the weak-coupling BCS limit 
to unitarity: a/ — *■ —00. In the unitarity limit there is 
a saturation of the effect of fermion interaction, prop- 
erly taken care of in the present set of equations, which 
seems to be appropriate to study the crossover [3J| from 
the weak coupling BCS limit to the molecular Bose limit. 

In the second part of the paper, by including the inter- 
action between a boson and fermion pair, we obtain a set 
of coupled equation for superfluid Bose-Fermi dynamics 
in the presence of a Bose-Fermi interaction. The model 
correctly describes the dynamics as the Bose-Bose repul- 
sion and Fermi-Fermi attraction are varied from weak- 
coupling to unitarity. This model is tested numerically 
as the different scattering lengths are varied over wide 
ranges covering the weak-coupling to unitarity transition 
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for a superfluid Bosc-Fcrmi mixture. 

In Sec. II we present the hydrodynamical equations for 
bosons and fermions and present the appropriate bulk 
chemical potentials valid in the weak-coupling limit as 
well as possessing the saturation in the strong-coupling 
unitarity limit consistent with the constraints of quantum 
mechanics. In Sec. Ill we derive mean-field equations for 
bosons and fermions consistent with the hydrodynam- 
ical equations for a large number of atoms. Next we 
derive the present model for interacting Bose-Fermi su- 
perfluid introducing a contact interaction between bosons 
and fermions. In Sec. IV we perform numerical calcula- 
tions for densities and chemical potentials of a trapped 
Bose-Fermi superfluid mixture to show the advantage of 
the present model. Finally, in Sec. V we present a brief 
summary and conclusion. 



II. SUPERFLUID HYDRODYNAMICS FOR 
BOSONS AND FERMIONS 



At zero temperature, for a large number of atoms, 
statical and dynamical collective properties of bosonic 
and fermionic superfluids are expected to be properly 
described by the hydrodynamical equations of superflu- 
ids [11, 13] ■ For bosons the hydrodynamical equations are 
given by [U, [3| 



d_ 

di 
d 



Ub + V ■ {ub-Vb) = 



at 



1 2 , ^ 

-mbVf^ +Ub + fJ,b[nb, ab) 



(1) 
0, (2) 



where Ubi^) is the external potential acting on bosons, 
•nib is the mass of a bosonic atom, nf,(r, t) is the local dcn- 
sityof bosons, and Vf,(r, t) is the local superfluid velocity 



sityot t 



The total number of bosons is given by 



Nb 



nb{r, t) 



dh 



(3) 



and the nonlinear term fib{nb,ab) is the bulk chemical 
potential of the bosonic system with Ub the Bose-Bose 
scattering length. Equation ([T]) is the equation of conti- 
nuity of hydrodynamic flow, while Eq. ^ is the equation 
of conservation of the linear momentum. Equation ([2]) es- 
tablishes the irrotational nature of the superfluid motion: 
V A Vb = 0, meaning that the velocity Vb can be written 
as the gradient of a scalar field. Equations ^ and ^ 
differ from the corresponding equations holding in the 
collisional regime of a non-superfluid system because of 
the irrotationality constraint. 

For superfluid fermions, the fundamental entity gov- 
erning the superfluid flow is the Cooper pair. In the 
case of a two-component (spin up and down) superfluid 
Fermi system one can introduce the local density of pairs 
rip(r,t) = nf{r,t)/2, where nf{r,t) is the total local 
density of fermions. The hydrodynamical equations of 



a fermionic superfluid 0, H, d, Q can be written as 
d 



dt 
d 



Up + V ■ [rip-Vp) = , 



impUp C/p + ^lp{np, af) 



(4) 

= , (5) 



where the nonlinear term fip(np,af) is the bulk chemi- 
cal potential of pairs with a f the attractive Fermi- Fermi 
scattering length. Here mp is the mass of a pair, that is 
twice the mass of a single fermion, i.e. rUp = 2m f. In 
addition, the trap potential Up{r) is twice the trap poten- 
tial Uf{r) acting on a single fermion, i.e. Up{r) = 2?7/(r). 
Note that the chemical potential ^p{np,af) is twice the 
total chemical potential ^f{nf,af) of the Fermi system, 
i.e. iJLp{np,af) = 2fif{2np,af). The total number Nf of 
fermions is 



TV 



(6) 



where Np is the number of pairs. 

It is important to stress that Eqs. ([T|) and ^ are 
formally similar to Eqs. Q and but the quantities 
involved are strongly different. In particular, the bulk 
chemical potential of bosons is completely different from 
the bulk chemical potential of Fermi pairs [l|, d, ■ 

In the full crossover from the small-gas-parameter 
regime to the large-gas-parameter regime we use the fol- 
lowing expression for the bulk chemical potential of the 
Bose superfluid 



I N 'i 2/3^/ 1/3 X 

IJ,b[nb,ab) = — n^' /(ji^' ab) , 
nib 



where 



fix) = 47r 



r5/2 



1 + 7^3/2 + /3a;5/2 



In the small-gas-parameter [x 0) regime, Eq. 
comes 



/(^) 



47r 



x + {a- 7)2;^/^ + 



(7) 

(8) 
be- 

(9) 



which is the analytical result for bulk chemical potential 
found by Lee, Yang, and Huang [3^, [s^l in this limit, 
provided that (a - 7) = 32/(3V7r). In Eq. dH) we shaU 
also choose /3 = ATTa/rj, with rj = 22.22; this will make 
the bulk chemical potential ([7]) satisfy the correct unitar- 
ity limit iibirib, ap) = 22.22h^Tv^^ jmb as established by 
Cowell et al. [3^. Hence Eqs. ([7]) with ([8]) can be made 
to satisfy the correct weak-coupling and unitarity limits. 
Next we need to determine the constants a and 7 consis- 
tent with (a — 7) = 32/(3-\/7r). We consider the following 
three choices for a and 7 from a possible many choices: 
a = 32j//(3V^), 7 = 32(i/- 1)/(3V7F) with (a) = 1.05, 
(b) v = 1.1, and (c) i' = 1.15. All three choices arc con- 
sistent with the unitarity and the Lee- Yang- Huang limits 
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The hydrodynamical equations are valid to describe 
equilibrium properties and dynamical properties of long- 
wavelengtli for both bosons and fermions. In particular, 
one can introduce P, a healing (or coherence) length 
such that the transport phenomena under investigation 
must be characterized by a length scale much larger than 
the healing length. As suggested by Combcscot, Kagan 
and Stringari [4J|, the healing length can be defined for 
bosons as 



FIG. 1: (Color online) The critical sound velocity v^^ of a 
uniform Bose superfluid as a function of Bose-Bose scattering 
length at and density nj. 



[36l . |37| . We shall present a critical numerical study of 
the three choices in the next section. 

In the full crossover from BCS regime (a/ — > —0) to 
unitarity (a/ ^ — c») we use the following expression for 
the bulk chemical potential of the Fermi superfluid [3§| 



-(67r2n,)V3g(2i/3„i/3„^) (ig) 



where 



1 



Sx 



1 



(11) 



where S and k are fitting parameters. The close agree- 
ment of the interpolation function g{x) for negative x 
values with the results of Monte Carlo calculation of en- 
ergy of a uniform gas of Fermi superfluid with the pa- 
rameters S = 47r/(37r2)2/3^ ^ = S/{1 - C), and ( = 0.44 
pol . I4H was established in Ref. [s^. The inclusion 
of the lowest-order term g{x) = 1 in Eq. ^TU\\ leads 
to the bulk chemical potential in the absence of Fermi- 
Fermi interaction (a/ =0) of a uniform gas. The next- 
order term g(x) = 1 + dx include known analytical re- 
sult in the small-gas-parameter regime as obtained by 
Lee and Yang [3a, \^ and by Galitskii [i^. The 
model (flO)) with ((TT|) provides a smooth interpolation 
between the bulk chemical potential in the BCS limit 
lip{np,af) — 2h^{6TT^np)^/^ /mp + lQirnpfi^af /rup [sl] for 

small rip^^af values and that in the unitarity limit [39j 

IJLp{np,af) = 2h^ idTr^nj,)^^^ C / iTT-v for large rip^^af values. 
Manini and Salasnich j2§| and Kim and Zubarev [2^ also 
proposed similar interpolation formulas on the basis of 
Monte Carlo data of a uniform Fermi superfluid. 

Here we deal with a Fermi gas in a spherical harmonic 
trap rather than a uniform Fermi gas. By a direct com- 
parison with the results of Monte Carlo calculation of 
a superfluid Fermi gas in a spherical harmonic trap, we 
shall show in Sec. Ill that Eqs. ((TO)) and (|lip still present 
a good approximation to the bulk chemical potential of 
the system, but now with a slightly different value of the 
parameters. 
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(12) 



where v^^ is the Landau critical velocity above which the 
system gives rise to energy dissipation. This critical ve- 
locity coincides with the first sound velocity and is given 
by 



irib dub 



(13) 



In Fig. [T] we plot the critical sound velocity vf of a 
uniform Bose superfluid for the present /i^ given by Eqs. 
([7]) and ([5]) with the parameter v = 1.1. The interesting 
feature of this plot is that w^'' = 0, when either rib or 
ab is zero. However, v^^ increases monotonically with 
Ub, whereas v^^ attains a constant value as ab increases. 
These features are present in the plot of Fig. [TJ However, 
the critical sound velocity v^^ calculated from /i^ given 
by Eqs. ([7]) and ([9]) increases monotonically with both 
Ub and ab- 

For superfluid fermions the healing length of Cooper 
pairs can be defined as 



n 



(14) 



where the critical velocity v'^ is related to the breaking 
of pairs through the formula 



(15) 



where |A| is the energy gap 0, S^]. We notice that in the 
deep BCS regime of weakly interacting attractive Fermi 
atoms (corresponding to |A| ^ ^p) Eq. approaches 
the exponentially small value v'^ = | A|/-\/mp/ip/2. 



III. NONLINEAR SCHRODINGER EQUATIONS 

The bosonic superfluid can be described by a CP [l].[25j 
complex order parameter 5'b(r, t) given by 



^b{v,t) = ^^^^)e'''^'''K (16) 
The probability current density i{v, t) is given by 



J = rib^b 



2imb 



(^^V^b-^-bV^-^). 



(17) 
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US Eq. (8); v=1.05 — 
US Eq. (8);v=l.l 
US Eq. (8);v=1.15 ■■■■ 

MGP Eq. (21) 

GP 

asymptotic: 22.22 
4 6 8 
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FIG. 2: (Color online) The Bose interpolation function f{x) 
given by Eq. ((8} with a = 32i//(30F), 7 = 32{u - 1)/(30F), 
and l3 = 47ra/22.22 with v = (a) 1.05, (b) 1.1, and (c) 1.15. 
The results of MGP [Eq. ^] and GP [f{x) = 4ttx] models 
together with the asymptotic value are also given. 



Equations (fTC|) and (fT7)) relate the phase 9b{r,t) to the 
superfluid velocity field v^, by the formula Q : 



V6 = — VOb . 

7716 



(18) 



The bosonic order parameter \E';,(r,t) is, apart from a 
normahzation [1, nothing but the condensate wave 
function 



Et{r,t) = (V'(r,t)) 



(19) 



that is the expectation value of the bosonic field operator 
The order parameter ^'b(r,<) satisfies the 
following generalized GP nonlinear Schrodinger equation 

i 



at 



2mf, 



+ [/fc + /-ifc(?l6, flb) 



(20) 



where the bulk chemical potential is given by Eq. ([7]). 
The use of f{x) = Attx in Eq. ([7]) leads to the GP equa- 
tion and the use of the two leading terms of Eq. ([9]), 
e-g-, 



fix) = 4tt 



32 



(21) 



in Eq. ([7]) leads to the modified GP (MGP) equation sug- 
gested by Fabrocini and Polls [IS] , which includes correc- 
tion to the GP equation for medium values of Bose-Bose 
scattering length Of,. The use of Eq. ([8]) in Eq. ([7]) leads 
to a nonlinear Schrodinger equation valid from the weak- 
coupling to the unitarity limit and is called the unitary 
Schrodinger (US) equation for bosons [2^ . 

In Fig. [2] we plot the various possibilities for the inter- 
polation function f{x); e.g., those corresponding to the 
US equation ^ with (a) 1/ = 1.05, (b) ly = 1.1, and (c) 
1/ = 1.15 [as suggested after Eq. ([§])], the MGP equation 



([2T|) . the GP equation f{x) = inx. The asymptotic limit 
lim2;^oo f{x) = 22.22 is also shown. All three choices 
(a), (b), and (c) represent smooth interpolation of f{x) 
between small and large x values. For small x this choice 
agrees with the MGP result ((2T|) and for large x with the 
asymptotic result. To test these three choices we actually 
solve the US Eq. , by imaginary time propagation af- 
ter a Crank-Nicholson discretization iZ] (detailed 
in the beginning of Sec. IV), for the trap parameters of a 
possible experimental set up for ^'^Rb atoms in a spher- 
ical trap and compare the results for energy with those 
obtained by Blume and Greene [i^ by diffusion Monte 
Carlo (DMC) method. Actually, we solved Eqs. (20) and 
(7) with rib normahzed to {Nb — 1) in place of Eq. (3), ap- 
propriate for a small number of bosons, cf. Eq. (2) of Rcf. 
[48] . The energy of the system is calculated through a nu- 
merically constructed energy functional from the present 
bulk chemical potential. The calculation is performed 
for a boson-boson scattering length a?, = 0.433^ where I 
is the harmonic oscillator length for various Nb and the 
results for energy are tabulated in Table [H From this 
table we find that the present US models always provide 
a much better approximation to the energies than the 
GP and MGP models. More results of energy for larger 
number of bosons may help in fixing the parameters of 
the present model more accurately. In the present paper 
we shall use the US model (b) with v = 1.1. 



TABLE L Ground state energies for different number Nt of 
bosonic atoms in a spherical trap and for ai, — 0.433 from 
a solution of GP equation, MGP equation and US equation 
for v = (a) 1.05, (b) 1.1 and (c) 1.15. The results for ener- 
gies obtained with two potentials of DMC calculation [4^ are 
quoted for comparison. Length and energies are expressed in 
oscillator units. 





DMC 


GP 


MGP 


US (a) 


US(b) 


US(c) 


3 


5.553(3); 5.552(2) 


5.329 


5.611 


5.570 


5.564 


5.558 


5 


10.577(2); 10.574(4) 


9.901 


10.772 


10.629 


10.608 


10.588 


10 


26.22(8); 26.20(8) 


23.61 


26.84 


26.24 


26.16 


26.07 


20 


66.9(4); 66.9(1) 


57.9 


68.5 


66.38 


66.07 


65.77 



The relationship between Eq. ([20]) and the hydrody- 
namical equations ([ij and ([2]) can be established by in- 
serting Eq. (fTB|l into Eq. (PD|) . After some straightfor- 
ward algebra equating the real and imaginary parts of 
both sides and taking into account Eq. (jl8|) . one finds 
two generalized hydrodynamical equations [l| 



d_ 

dt 
d 



rib 



V • (rifcVb) = 



at 



V 



rubvl 



2mb 
fJ-birib, ab) 



Ub 

0, 



Ub 



(22) 



(23) 
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FIG. 3: (Color online) The Fermi interpolation function 
g{x) given by Eq. ([TTJ with (i) 5 = 47^/(37^^)^''^ and (u) 
207r/(37r^)^/^ and k = 5/0.56. The MDF results given by Eq. 
(|29|) and the asymptotic value are also shown. 



(24) 



which include the quantum pressure term 



T, 



QP 



2mf, 



which depends explicitly on the reduced Planck constant 
h. Neglecting the quantum pressure term, one gets from 
Eqs. ((22|) and ((23|) the classical hydrodynamical equa- 
tions (III) and Equation is the continuity equa- 
tion whereas Eq. establishes the irrotational nature 
of hydrodynamic flow. 

Similarly, the fermionic superfluid can be described by 
a DF [sii complex order parameter 'I'p(r, t) of pairs, with 
a modulus and a phase Op{r, t) such that 



4'p(r,<) = ^np{v,t) 



(25) 



The expression for the probability current density of 
fcrmions leads to the following expression for the velocity 
field Vp of pairs [1| : 



(26) 



Note that the mass of a pair TOp, and not that of a 
fcrmion, appears in the denominator of the phase- velocity 
relation. The fundamental entity responsible for super- 
fluid flow has the mass mp of a pair of fermions. For 
paired fermions in the superfluid state, the order param- 
eter is, apart from a normalization [1, 0, S^, the conden- 
sate wave function of the center of mass of the Cooper 
pairs: 



Hp(r,t) = (V'T(r,t)i/'i(r,t)) 



(27) 



that is the average of pair operators, with ip^lr^t) the 
fermionic field operator with spin component a =t, i 0j 
[23} . The order parameter 4'p(r,t) satisfies the following 
zero-temperature nonlinear DF Schrodinger equation |23l . 



9 T 



2m„ 



Up + fip{np,af) 



(28) 
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FIG. 4: (Color online) (a) Energy of a superfluid Fermi gas 
in a spherical trap in oscillator units vs. Fermi- Fermi scatter- 
ing length af also in oscillator units obtained from the solu- 
tion of Eq. (ggi for the choices S = (i) 47r/(37r^)^/^ and (ii) 
207r/(37r^)^/^ and from the fixed node Monte Carlo (FNMC) 
calculation [s^. Iss!]. (b) Energy of a superfluid Fermi gas in 
a spherical trap in oscillator units vs. number of atoms Nf 
in the unitarity limit af — > —00 obtained from a solution of 
Eq. ((28ll, from FNM C tsoll and Green function Monte Carlo 
(GFMC) calculations [5J| , and from local density approxima- 
tion (LDA). 



where the bulk chemical potential /ip(np, ay) is given by 
Eq. dTO]). The use of .g(a;) = 1 in ^ip{np,af) of Eq. ^ 
leads to the lowest order DF equation. If we use the next 
order solution of Eq. (fTTjl in Eq. (fTO|) . e. g., 



g{x) = 1 + 5x, 



(29) 



we obtain a modified density-functional (MDF) equation 
with corrections for medium values Fermi-Fermi scatter- 
ing length a / . The use of Eq. (fTTjl in Eq. (fTO|) leads to a 
DF equation valid in the weak-coupling BCS to unitarity 
limit and the resultant nonlinear equation (|28p will be 
called the unitary Schrodinger (US) equation for fermion 
pairs. 

In Fig. [3] we plot the various possibilities for the inter- 
polation function g{x)] e. g., that corresponding to Eq. 
(|TT|). the MDF equation ([29)) for two different values of 
5: (i) 47r/(37r2)2/3 and (ii) 207r/(37r2)2/3. The parameter 
K is always taken as k = (5/(1 - C), C = 0-44 [13, EH. 
The asymptotic limit lima;^oo g{x) = 0.44 is also shown 
in Fig. [3l The expression ([TT|) represents a good approx- 
imation of g{x) for small and large |a;|. For small |a;| this 
choice agrees with the MDF result (|29p and for large \x\ 
with the asymptotic result for both choices of 6. 
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To test the above two choices (i) and (ii) of the param- 
eter S in Eq. (jlip for a Fermi superfluid in a spherical 
trap, wc calculate the energy of the system for different 
values of the scattering length a f and number of atoms 
Nf by solving directly Eq. ((28|) by imaginary time propa- 
gation after a Crank-Nicholson discretization [i^ 0, [irj 
(detailed in the beginning of Sec. IV). We also compare 
the results with those obtained by the fixed-node Monte 
Carlo calculation (FNMC) [13, [5l|, [H, HI] . The results 
of our investigation are shown in Figs. SI In Fig. 3] 
(a) we plot energy E vs. |a/| for Nf = 4 and 8 ob- 
tained from a solution the US equation with the 
choices (i) 47r/(37r2)2/3 and_(ii) 207r/(37r2)2/3 for 6 and 
the FNMC data of Ref. [52|, \5^. The present hydro- 
dynamical formulation is expected to be good for for a 
large number of fcrmions [see, Fig. 0] (b)], yet for a small 
number Nf = 4 the result is reasonable. In Fig. [4] (b) 



we plot energy E/N^^ vs. Nf'"^ at unitarity a/ —> — oo 
obtained from a solution of the US equation (pS)) with 
choice (ii) 207r/(37r^)^/^ for 6. In this limit both choices 
of S lead to the same energy. The results of FNMC [50| 
and Green function Monte Carlo (GFMC) [11] calcula- 
tions as well as of LDA are also shown in Fig. 2] (b). 
The LDA result is analytically known in this case as 
E{N) = (3iV/)4/37J/4,^ = 0.44 [33 

between these two variables explicit in the LDA result. 
From Figs. 5] we find that the results with choice (ii) 
207r/(37r2)2/^ of S agree well with the Monte Carlo data 
in both cases and this choice will be used in the present 
study. The Monte Carlo data clearly favors the present 
model over the LDA. 

Again, inserting Eq. into Eq. after some 

straightforward algebra equating the real and imaginary 
parts and taking into account Eq. (QH) , we find two gen- 
eralized hydrodynamical equations for the Fermi super- 
fluid 



2/3 



We plot E/N^^ 



vs. Nf' in Fig. [11(b) because of the linear correlation 



d 



V • (jipVp) = , 
V 



2m.n 



+ t^p{np,af) 



= 



which include the quantum pressure term 



rQP 

p 



2toe 



(30) 



(31) 



(32) 



involving h. Neglecting this term leads to the classical hy- 
drodynamical equations ^ and ([5]) of superfluid Fermi 
pairs. This specific quantum pressure term is a conse- 
quence of the proper phase- velocity relation ((26)) with 



the factor irip in the denominator [0]. Had we taken a 
different mass factor in the denominator, e.g. m/, a dif- 
ferent quantum pressure term would have emerged. The 



present choice is physically motivated by Galilei invari- 
ance [3, [H, [3l| and leads to energies of trapped Fermi 
superfluid in close agreement [sl] with those obtained 
[sol . [5^ from the Monte Carlo methods. 

We stress that, from the point of view of DF theory, 
the quantum pressure terms for superfluid bosons and 
fcrmions correspond to gradient correction (also called 
surface corrections) to the local density approximation 
(LDA), where LDA gives exactly the classical hydrody- 
namical equations obtained by setting gradient kinetic- 
energy term to zero. The gradient term, which can also 
be seen as a next-to-leading contribution in a momen- 
tum expansion to the classical superfluid hydrodynamics 
[5^, takes into account corrections to the kinetic energy 
due to spatial variations in the density of the system, 
and it is essential to have regular behavior at the sur- 
face of the cloud where the density goes to zero [STj . 
In the case of superfluid bosons, the quantum-pressure 
term, Eq. (|24L is exactly the one of the familiar cubic 
CP equationfil, . In the case of superfluid fermions 
the gradient term is consistent with the motion of paired 
fermions of mass rrip = 2to/. For normal fermions var- 
ious authors have proposed and used different gradient 
terms [13, HI, [H, [13, |6l| . For superfluid fermions we are 
using the familiar von Weizsacker term [1^ in the case of 
pairs. In fact, in our approach, the explicit form of the 
quantum-pressure term, Eq. (|32p . is strictly determined 
by the Galilei- invariant relation, Eq. (|26p . between phase 
and superfluid velocity [2, [H, [3l|. Finally, we observe 
that Eq. ([5^ is practically the same as that recently ob- 
tained in Ref. [B^] in the case of a superfluid Fermi gas 
at unitarity using an e expansion around c? = 4 — e spatial 
dimensions. 

The generalized hydrodynamical equations (j22p and 
for bosons and (|30p and ([?T|) for fermions are 
thus completely equivalent to the respective nonlinear 
Schrodingcr equations ^U\i and 
velocity flelds, Vf,,Vp 0. Eqs. 
stationary versions of Eqs. (PO)) 
superfluid densities: 



In the limit of zero 
(HSl) and dni) yield the 
and (pS)) describing the 



2m, 



■V^ + Ub+ fibin-b, at) 



2m, 



-V^ + U.p+ flpiUp, Qf) 



rib , (33) 



(34) 



showing the agreement between hydrodynamical and 
nonlinear equations for superfluid bosons and fermions. 
Here fXbo and are the chemical potentials for bosons 
and fermions in the stationary state. 

Although Eq. (j34[) is written in terms of Fermi pair 
density Up (and we use this equation in the following), 
an equivalent equation for Fermi density n/ (=2np) can 
be written as 



8m t 



V'+ Uf+ iif{nf,af) 



Uf , (35) 



tJ-finf^af) 



2m i 



{iir^nff'^g{n)'\f) , (36) 
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where we have used /ipo = 2^yo, Up = '^Uf, and 
IJ,p{np,af) = 2fj,f{nf,af). 

Equations pU)) and (|28|) are the Eulcr-Lagrange equa- 
tions of the following Lagrangian density 



ih 



d 



d 



C- — Vl/* \T/ ■ — Vl/ Vl/ 



2m 1 



(37) 



with j = &, p, (here ap = aj, the Fermi- Fermi scattering 
length) where 



(38) 



is the bulk energy per particle of the superfluid. 

Finally, we introduce a Lagrangian density for Bose- 
Fermi interaction to ([57)) of the form 



Cbp = GhJ'fhl^l'^.W 



'"bp I 



(39) 



where Gbp = inh^abf /nibf, and mt,f is the Bose-Fermi 
reduced mass and Ubf the Bose-Fermi scattering length. 

We note that there are no coupling terms between the 
hydrodynamical equations (|23|) and PT|) involving the ve- 
locity of two superfluid components. The only coupling 
in the Lagrangian density (|39p involves the density of 
the two superfluid components. However, the present hy- 
drodynamical equations and (|3ip are a special case 
of the three-fluid (two superfluids and a normal fluid) 
hydrodynamical formulation of Landau and Khalatnikov 
[631 (appropriate for the ^He-^He mixture). The Galilei- 
invariant Landau-Khalatnikov formulation in its general 
form also has a coupling term involving two superfluid ve- 
locities (not considered in this paper), as discussed (G^I 
by Andreev and Bashkin and also by Ho and Shenoy. 
However, the present static model defined by Eqs. ([55]) 
and (pi]) is obtained by setting all velocities equal to zero, 
and hence, these terms depending on superfluid velocities 
would leave the present model unchanged. 

In this paper we consider ab / to have not too large pos- 
itive values only (moderately repulsive Bose-Fermi inter- 
action.) There are other domains of Bose-Fermi interac- 
tion which require special attention. If this interaction 
is attractive and strong enough, then the ground state 
is very different, e.g., the composite fermions f = bf 
are created at higher temperatures and at zero temper- 
ature we have Cooper pairs (or local pairs) of the type 
< / / > consisting of two elementary fermions and two 
elementary bosons . 

The complete Lagrangian density given by Eqs. (j37p 
and leads to the following set of coupled equations 
for the Bose-Fermi-Fermi mixture made of superfluid 
bosons and two-component superfluid fermions 



th—'^b = 

at 



2mfc 



+ Ub + ^ib{nb, ab) + Gbp\'^pY 



*6, 

(40) 
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FIG. 5: (Color online) The normalized densities nj{r) = 
with j = b,p from a solution of Eqs. (|40|l and (|4H) for the 
US, MGPDF, and GPDF models with Nb = 10000, A^p = 
10000, a;,/ = 100 nm, aj — 0, and (a) at = 200 nm, and (b) 
at = 2000 nm. In this plot / n■j{r)d^^r = 1. 



^n-^p = 



2m„ 



+ Up + fJ-piup, af) + Gfcpl^'hl^ 



p ' 
(41) 



valid from weak-coupling to unitarity for both bosons and 
fermions. The coupled set of equations and ([^T|) 
with bulk chemical potentials /i(,(nfc,af,) and ^p(np,af) 



given by Eqs. ^ and (O (with i/ = 1.1) and Eqs. ([TI 
and (HI]) [with 5 = 207r/(37r2)2/3^K = ,5/(1-0, C = 0.44], 
respectively, is the US equation for the interacting super- 
fluid Bose-Fermi mixture and is the required set. How- 
ever, when the bulk chemical potentials yUb(nb,at,) and 
/ip(np,a/) are approximated using Eqs. ((2T|) and (|29p . 
the set of equations and (|¥T|) will be termed modi- 
fied Gross-Pitaevskii density-functional (MGPDF) equa- 
tion including some correction for medium values of scat- 
tering lengths flf, and af. Finally, if we approximate 
f{x) = Attx and g{x) = 1 in Eqs. (gO]) and (gT]) with 
Eqs. ^ and (fTT|) . the resultant model is termed Gross- 
Pitaevskii density-functional (GPDF) model. We present 
a comparative numerical study of these three models — 
US, MGPDF. and GPDF - in the next section. 



IV. NUMERICAL RESULTS 

Here we solve numerically the coupled set of US Eqs. 
(|40)) and (|4T|) with bulk chemical potentials fib and fj,p 
given by Eqs. ([5]) and (fTTj) and compare the results so 
obtained with those of the MGPDF and GPDF mod- 
els. We numerically solve the US, MGPDF, and GPDF 
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FIG. 6: (Color online) The normalized densities nj{r) = 
with j = b,p from a solution of Eqs. (|40|l and (|41|l for the 
US, MGPDF, and GPDF models with Nt = 100000, Np = 
100,06/ = 100 nm, at — 5 nm, and (a) a/ = —100 nm, and 
(b) atf = —300 nm. In this plot J nj{r)d'^r — 1. 



partial differential equations by discretizing them by the 
semi-implicit Crank-Nicholson algorithm with imaginary 
time propagation [i^, [4^, \^ . The space and time steps 
used in discretization were typically 0.05 and 0.001 re- 
spectively. 

Although we are not interested here in fitting experi- 
mental data, we perform numerical calculation for vari- 
ables for a possible experimental set up [3 of *^Rb- 
^"K mixture in spherical harmonic traps satisfying Ub = 
mijLolr^ /2 = ■mfUj'jr'^/2 = Uf with an oscillator length 

I = y^h/ (mfoWb) = 1 /im for ®^Rb atoms. The oscillator 
length for is modified according to the trap. The 
fermion pair trap in Eq. (j4ip is Up = nifUj'jr'^. 

In Figs. O (a) and (b) we plot densities for the super- 
fluid Bose- Fermi mixture as calculated from Eqs. ((40l) 
and gl]) for Nt = 10000, iVp = 10000, a/ = 0,06/ = 100 
nm with (a) ai, = 200 nm and (b) ai, = 2000 nm for 
the US, MGPDF, and GPDF models. For smaU at both 
MGPDF, and GPDF models are good approximations to 
the bosonic density of the US model (not shown in fig- 
ure). But as Qb increases, the MGPDF model is a better 
approximation to the bosonic density of the US model 
than the GPDF model as can be seen from Fig. [5] (a). 
However, for very large Ub, neither the MGPDF nor the 
GPDF model provides the saturation of bulk chemical 
potential for bosons as in the US model. Consequently, 
the MGPDF, and GPDF models acquire large bosonic 
nonlinearity compared to the US model and hence pro- 
duce bosonic density extending to larger region in space. 
As Oft is increased, the MGPDF model yields rapidly di- 



FIG. 7: (Color online) (a) Chemical potential fito vs. at of 
Eqs. ^ and ^ for the US, MGPDF, and GPDF models 
for Nt ^ Np ^ 10000, a/ = 0, at/ = 100 nm and (b) chemical 
potential ^^po vs. a/ for the US, MGPDF, and GPDF models 
for Nb = 100000, Np = 100, at = 5 nm, flb/ = 100 nm, respec- 
tively. The chemical potentials are expressed in oscillator unit 



vergcnt bosonic nonlinearity (compared to the GPDF 
model) and hence produces bosonic density inferior to 
the GPDF model as can be seen in Fig. [5] (b). 

In Figs. [6] (a) and (b) we plot superfluid densities 
for Nb = 100000, iVp = 100, Ob/ = 100 nm, = 5 nm 
with (a) af = —100 nm and (b) cif = —300 nm for the 
US, MGPDF, and GPDF models. For small |a/| the 
fermionic densities of both MGPDF, and GPDF models 
are good approximation to the US model. But as |a/| 
increases, the GPDF model continues with BCS nonlin- 
earity for af = 0, whereas the fermionic nonlinearity of 
the MGPDF model acquires excessive attraction. This is 
clear from Figs. [S](a) and (b) showing a peaked fermionic 
density for the MGPDF model for large |a/| correspond- 
ing to increased attraction. For very large |a/|, the MG- 
PDF model does not provide the saturation of fermionic 
bulk chemical potential as in the US model. We note 
that compared to the study of Figs. [5] with comparable 
number of bosons and fermions, in the study of Figs. [6] 
we have much reduced number of fermions compared to 
bosons and the system has undergone a mixing-demixing 
transition (20l . [6a | expelling the fermions from the central 
to the peripheral region. 

Next we calculate the chemical potentials fibo ^^nd 
of the stationary superfluid Bose-Fermi mixture as de- 
scribed by Eqs. §^ and In Figs. [7] (a) and (b) we 
plot these chemical potentials under the situations de- 
scribed in Figs. 0] and [SI e.g., in the (a) weak-coupling 
BCS limit when varies from a small to a large value 



9 



and (b) weak-coupling GP limit when |a/| varies from 
a small to a large value, respectively. In the first case 
the fermionic chemical potential fipo remains practically 
constant over the full range of variation of and we plot 
only the bosonic chemical potential /if,o in Fig. [7] (a). In 
the second case bosonic chemical potential Hbo remains 
practically constant over the full range of variation of ay 
and we plot only the fermionic chemical potential fipQ. 
From Fig. [7] (a) we find that, although for small val- 
ues of Of,, the MGPDF model produces results for ^i,o in 
close agreement with the US model, for medium values 
of ab the GPDF model produces results closer to the US 
model. However, for larger values of at the US model 
should show saturation, whereas the GPDF chemical po- 
tential fibo should increase monotonically with ab and 
should be larger than the the chemical potential of the 
US model (not shown in figure). From Fig. [21(b) wc find 
that none of the MGPDF or GPDF models produces rea- 
sonable result for /i^g except for very small values of \af\. 
The GPDF model does not have a dependence on ay and 
hence produces a constant result for fipQ. 

V. CONCLUSION 

The dynamics of cold trapped atoms (both bosons and 
fermions) should be correctly handled by a proper treat- 
ment of many-body dynamics or through a field-theoretic 
formulation. However, both approaches could be much 
too complicated to have advantage in phenomcnologi- 
cal application. This is why mean- field approaches are 
widely used for phcnomcnological application. In the 
weak-coupling limit, for bosonic atoms the dynamics is 
handled through the mean-field GP equation [l| , whereas 
for fermionic atoms the LDA (local density approxima- 
tion) formulation is widely used The LDA approach 
(often called the Thomas- Fermi approximation) lacks the 
kinetic energy gradient term (quantum pressure term). 
Here we introduce a quantum pressure term consistent 
with the correct phase-velocity relation of the superfluid 
fermions assuring the Galilei-invariance of the model [1^ . 
This modified LDA equation for fermions and the GP 
equation for bosons, both valid in the weak-coupling 
limit, are both shown to be equivalent to the hydro- 
dynamical equation for superfluid flow. In the strong- 



coupling unitarity limit as the bosonic and fermionic 
scattering lengths ab and |a/| tend to infinity, both the 
bosonic and fermionic interactions are to exhibit satura- 
tion due to constraints of quantum mechanics. We intro- 
duce proper unitarity saturation effects in these equa- 
tions. In addition in our Bosc- Fermi formulation we 
introduce a not too strong contact interaction between 
bosons and fermions. The resultant superfluid Bose- 
Fcrmi dynamics is described by a set of coupled partial 
differential equation, which is the principal result of this 
paper. In the time-dependent form with a Bose-Fermi 
interaction this set is described by Eqs. (HO]) and (pij) 
with the stationary version without Bose-Fcrmi interac- 
tion given by Eqs. ^ and ((M)) . 

The present model for the superfluid Bose-Fermi mix- 
ture is termed the US model, valid in both weak- and 
strong-coupling unitarity limits for bosons and fermions. 
We also consider a model called the GPDF model where 
the boson dynamics is handled by the GP Lagrangian and 
the fermion dynamics is handled by the LDA formulation 
with a proper kinetic energy gradient term. Finally, we 
also consider another model called the MGPDF model in- 
cluding lowest order correction to the GPDF model due 
to the non-zero bosonic and fermionic scattering lengths 
Ub and \af\. However, the MGPDF model does not pro- 
vide a saturation of the nonlinear interaction in the uni- 
tarity limit when ab and |ay| tend to infinity. In our 
numerical studies of a ^^Rb-'*'^K mixture in a spherical 
harmonic trap, we find that for medium interactions the 
MGPDF model produces improved results. However, for 
stronger interactions the full US model should be used. 
Although, slightly complicated for analytic applications, 
the present US model is no more complicated than the 
usual GP and LDA equations for numerical treatment 
and is expected to find wide applications in future stud- 
ies of Bose-Fermi superfluid mixture. 
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